Phase transitions in crowd dynamics of resource allocation 
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We define and study a class of resources allocation processes where gN agents, by repeatedly 
visiting N resources, try to converge to optimal configuration where each resource is occupied by 
at most one agent. The process exhibits a phase transition, as the density g of agents grows, from 
an absorbing to an active phase. In the latter, even if the number of resources is in principle 
enough for all agents (g < 1), the system never settles to a frozen configuration. We recast these 
processes in terms of zero-range interacting particles, studying analytically the mean field dynamics 
and investigating numerically the phase transition in finite dimensions. We find a good agreement 
with the critical exponents of the stochastic fixed-energy sandpile. The lack of coordination in the 
active phase also leads to a non-trivial faster- is-slower effect. 

PACS numbers: 05.70.Fh, 89.65.-s, 87.23. Ge 



I. INTRODUCTION 

A general question that naturally arises in the study of 
social systems is how collective behaviors can emerge 
from the interactions among individuals. In spite of 
the inherent complexity of these phenomena, simplified 
mathematical models that assume simple automatic re- 
sponses of individuals to stimuli can reproduce non- 
trivial effects in the observed behavior [1]. Statistical 
mechanics had a certain success in the use of coarse- 
grained models of physical systems to connect the micro- 
scopic dynamics with the macroscopic behavior, and its 
techniques and concepts are starting to be fruitfully ap- 
plied to understand social dynamics 2-5]. In particular, 
the phenomena observed in the crowd dynamics, from 
pedestrians flows @ to vehicular traffic pj, have been 
recently subject to quantitative measures. Observations 
range from bottleneck oscillations, lanes and stripes for- 
mation, intermittent flows, waves, turbulence [81 , faster- 
is-slower and freezing- by- heating effects [l(|. Inter- 
estingly, this extremely broad class of collective phenom- 
ena is explained at semi-quantitative level with the use 
of models of interacting particles or granular fluids. This 
shows that even if interactions among individuals in these 
settings are rather simplified and mechanical, they still 
lead to a certain level of collective coordination. 

In this paper we show how a crowd dynamics in a re- 
sources allocation context can give rise to a phase tran- 
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sitions between an active and an absorbing phase. The 
class of models we discuss is inspired by the Kolkata Paise 
Restaurant (KPR) problem [Ill4l4| . which we generalize 
and recast in the broader context of zero-range interact- 
ing particle models fl5j . 

The KPR is a repeated game played by a large num- 
ber (gN, g being a real number) of agents with equally 
shared past information. Every evening they try to get 
the best service in one of the N restaurants of commonly 
agreed ranks, each providing food for one person only. It 
serves as a paradigm for a problem of resource utiliza- 
tion, where agents learn from their actions to maximize 
the effective utilization of available resources. At unit 
density (g — 1), a simple random choice algorithm (Tl| 
leads (in unit convergence time) to a occupation of a 
fraction of / = 1 — er 1 ~ 0.63 only of the resources, 
which of course falls much short of a fully efficient us- 
age of them. As seen in various earlier studies [pj - flij . 
apparently smarter strategies fail to yield better results. 
However, a stochastic strategy .14] that maintains a naive 
tendency to stick to an agent's past choice, with proba- 
bility decreasing with the crowd size, leads to an efficient 
utilization fraction / of about 0.796 [l4| and it converges 
to the above within a time (iteration) independent of N. 
An extension of this strategy to the Minority game [l6[ 
problem (with knowledge of both the sign and size of 
the population difference among two alternative choices), 
gives the optimal solution (with fluctuation bounded by 
log N and convergence time bounded by log(log N) [171]). 
Our focus will be on the nature of the collective behavior 
rather than on the efficiency of individual behavior. For 
the latter, we refer to applications of game theory and 
adaptive learning, as e.g. in |18l - l22| . We note, in pass- 
ing, that the relation between the degree of individual 



2 



rationality and the efficiency of collective allocation can 
be a non-trivial one, as shown e.g. in Refs. [23l |24|| . 

Here we will study the general problem with N restau- 
rants and gN agents, where g, the density, is a fixed ex- 
ternal parameter of the dynamics. Recasting the problem 
in terms of zero-range interacting particles, we report a 
phase transition from a frozen phase with satisfied agents 
to an active phase with unsatisfied ones at a critical den- 
sity g c and perform extensive numerical simulations as 
well as some analytical calculations to understand its fea- 
tures, finding a good agreement with the exponents of 
stochastic fixed-energy sandpiles [25]. The study of the 
relaxation properties of the frozen phase reveals an in- 
teresting faster- is-slower effect: This consists in the gen- 
eral observation that, from vehicular traffic to logistic, 
a higher level of coordination can arise from strategies 
which involve a slower dynamics (e.g. waiting), that how- 
ever speed up performances. 

The paper is organized as follows: In Sec. UH we intro- 
duce the class of models under study. Next, in Sec. IIIII 
we present the results of numerical simulations and in 
Sec. IIVI we provide the analytic treatment for the mean 
field case. We conclude with a summary and some dis- 
cussions in Sec. |V] 



II. THE MODELS 

Inspired by the KPR problem [TTJ, [l4| with N individuals 
competing for N restaurants (each serving food to one 
person each), we propose a more generic and fundamental 
stochastic occupation problem with exclusion. 

The rank ordering amongthe restaurants, as in the 
original KPR problem [TJ (3] is dropped here and we 
consider in general gN agents, where g, the density, is an 
external fixed parameter of the dynamics. We will refer 
from now on to individuals as particles and restaurants 
as the sites or nodes of an underlying graph. In these 
terms the original problem is defined in a fully connected 
geometry. 

A particle moves from the site i it occupies to a ran- 
domly chosen neighbor j with a rate that depends only on 
the number of particles that is present on it v(rii). This 
is by definition a zero range process (see [l5| for a review) 
which allows us to conclude that the stationary probabil- 
ity distribution of the number of particles per site factor- 
izes in terms of single site functions. Given the nature of 
the problem, we will consider models with v(l) = 0, i.e. 
costumers are happy while alone, and v(n+l) > v(n), i.e. 
repelling particles. Given this expression for the rates, at 
low densities (<? < 1) there are dynamically frozen con- 
figurations with sites filled by single particles. On the 
other hand, for high densities (g > 1) a finite fraction of 
sites - that we shall call 'active' - will have multiple oc- 
cupancy. We will show that there is a transition between 
these two phases that occurs at a certain density g c < 1 . 
Notice that, in principle, the process is ergodic so every 
configuration of is accessible. Therefore, for g < 1, the 



process will sooner or later visit an absorbing configura- 
tion where n.; < 1 for all sites i. Still, when N is large 
enough and g > g c , frozen configurations will be visited 
extremely rarely, over time-scales which are much longer 
than those which are accessible to numerical simulations. 
We choose as the order parameter the steady state den- 
sity of active sites p a (density of sites having n > 1). 
So the absorbing phase corresponds to p a = whereas, 
above some density g c the steady state shows a persistent 
dynamics with a non-zero value of the order parameter 

(Pa > 0). 

We will analyze in particular two models: 

(A) v(n) =1-1 

(B) v{n) = (l-p)9(n - 1). 

We will refer from now on to a parallel dynamics in 
which a simultaneous update of the nodes is done at each 
time step, i.e. agents' actions are simultaneous akin to a 
repeated game setting. 

It will be pointed out later that a sequential update 
in which at each step a randomly chosen particle jumps 
with some probability leads here to the trivial result that 
g c = 1. This corresponds to the limit p — > 1 of model B. 

The model A implements the stochastic crowd avoiding 
strategy of the original KPR problem [l4j|. Here, if the 
site k has rifc > 1 particles, each of them stay back with 
probability 1/nfe in the next time step, otherwise it jumps 
to any of the neighboring sites. 

In the model B an external parameter p is introduced 
that represents the "patience" of costumers to over- 
crowded conditions. Here, if the site k has rik > 1 par- 
ticles, each of them stays with probability p in the next 
time step, otherwise it jumps to any of the neighboring 
sites. 

The model B is in practice a kind of fixed energy sand- 
pile [25| , but the study of its dynamics as a function of 
the parameter p will reveal an interesting faster-is slower 
effect related to the relaxation time of the frozen phase. 

We finally point out that a waiting choice can be ratio- 
nal from the point of view of game theory: the agents in 
overcrowded sites could wait simply because they expect 
that others are leaving them alone. 

In the next section we report a detailed analysis of the 
phase transition of both models with the use of numerical 
simulations in a fully connected geometry, a Id chain 
and a 2d square lattice. Then, we perform analytical 
calculations for the mean field case. 



III. RESULTS FROM NUMERICAL 
SIMULATIONS 

We measure the times required to reach the steady state 
below and above g c . Below g Cl the order parameter p a 
reaches a value p a = in the steady state. Above g Cl 
order parameter p a evolves to a stationary state and fluc- 
tuates around a mean value p° a (> 0). The system has 
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FIG. 1. Simulation results for mean field case, g c = 0.7502 ± 0.0002. (a) Variation of steady state density p a of active sites 
versus g — g c , fitting to /3 — 0.98 ± 0.02. The inset shows the variation of p a with density g. (b) relaxation to absorbing 
state near critical point for different system sizes, the inset showing the scaling collapse giving estimates of critical exponents 
a — 1.00 ± 0.01 and z' = 0.50 ± 0.01. (c) Scaling collapse of p a (t). The inset shows the variation of p a (t) versus time t for 
different densities g. The estimated critical exponent is vu = 1.00 ± 0.01. The system sizes N are mentioned. 




FIG. 2. Simulation results for 2-d case, g c = 0.8827±0.0002. (a) Variation of steady state density p a of active sites versus g — g c , 
fitting to P = 0.68 ± 0.01. The inset shows the variation of p a with density g. (b) relaxation to absorbing state near critical 
point for different system sizes, the inset showing the scaling collapse giving estimates of critical exponents a — 0.42 ± 0.01 and 
z = 1.65 ± 0.02. (c) Scaling collapse of p a (t). The inset shows the variation of p a (t) versus time t for different densities g. The 
estimated critical exponent is i/n = 1.24 ± 0.01. The simulations are done for square lattices of linear size L (N = L 2 ). 



persistent dynamics in this phase. The evolution of the 
order parameter is exponential away from g Cl and can be 
expressed as 



Pa(t) 



Pa 



for g > g c , and 



Pa(t) 



1 - e 



-t/T 



-t/r 



(1) 



(2) 



for g < g Cl where r is the relaxation time. We will de- 
note the asymptotic value of the order parameter as p a 
hereafter. Near the critical point (g — g c — > 0+), we find 
Pa ~ (g — 9c)^ where j3 is the order parameter exponent, 
and t ~ (g — g c )~ v ^ . Generally p a (t) obeys a scaling form 



Pa(t) 



I \ '- 



{g-g c )-^~L*, (3) 



where a and z are dynamic exponents and L denotes 
size of the system. We then get (3 = v\\a by comparing 
Eq. dTJ, Eq. (0) and Eq. (|3|) when t/r is a constant for 
t — >• oo . We study numerically the time variation of p a (t) 
and measure the exponents by fitting to the above scaling 
relation. 



A. Model A 

1. Mean Field case 

For the mean field case, we have studied systems of N — 
10 6 sites, averaging over 10 3 initial conditions. We get 
g c = 0.7502±0002. The scaling fits of p a (t) for different g 
values (see Fig.Q} give /3 = 0.98±0.02, z' = 0.50±0.01 (if 
we assume N = L 4 and using Eqn. ([3]), we get a relation 
z = Az' and therefore z = 2.0 ± 0.04), V\\ = 1.00 ± 0.01, 
a = 1.00±0.01. It may be noted that these independently 
estimated exponent values satisfy the scaling relation /3 — 
v\\a well. 



2. Lattice cases 

We studied the same dynamics in 1-d and 2-d. For a lin- 
ear chain in 1-d, we took N = L = 10 4 and averaged over 
10 3 initial conditions. For 2-d we consider square lattice 
(JV = L 2 ) with L — 1000 and averaging over 10 3 initial 
conditions. Periodic boundary condition were employed 
in both cases. 
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FIG. 3. Phase diagram for the generalized model in the (g,p) 
plane, showing the phase boundaries separating the active 
and absorbing phases in 1-d, 2-d and mean field cases. The 
active phases are on the right of the phase boundaries while 
the absorbing phases are on the left in the respective cases. 
The system sizes are TV = 10 s for mean field, 1000 x 1000 for 
2-d, and 10 4 for 1-d. 



2. Lattice cases 

We studied the same dynamics in 1-d and 2-d. For a 
linear chain in 1-d, here also we took N = L = 10 4 
and average over 10 3 initial condition. For 2-d we took 
1000 x 1000 square lattice with L = 1000 and averaging 
over 10 3 initial conditions. 

(a) For 1-d, for the case p — 0, we find g c — 0.89 ± 0.01, 
with (3 = 0.42±0.01, z = 1.55±0.02, ^ = 1.90±0.02 
and a — 0.16 ±0.01 (Fig. @|. The phase boundary in 
(g,p) is nonlinear: it starts from g c = 0.89 ± 0.01 at 
p = (Fig. gj to p = 0.43 ± 0.03 at g = 1 (Fig. E}. 
Thus, one can independently define a model at unit 
density (g = 1) and calculate the critical probability 
p c for which the system goes from an active to an 
absorbing phase. 

(b) For 2-d, for the casep = 0, we find g c = 0.683±0.002, 
with (3 = 0.67±0.02, z = 1.55±0.02, v\\ = 1.20±0.03 
and a — 0.42 ± 0.01. The phase boundary seems 
nonlinear, from g c = 0.683 ± 0.002 for p = (Fig. \5§ 
extending to g c = 1 at p = 1. 



(a) The model is defined in 1-d as follows: The parti- 
cles are allowed to hop only to their nearest neighbor- 
sites, and each particle can choose either left or right 
neighbor randomly. We find g c = 1 and hence the 
phase transition is not very interesting. 

(b) In the 2-d version of the model, we consider square 
lattices and the particles are allowed to choose one of 
the 4 nearest neighbors randomly. For N = 1000 x 
1000, we get g c = 0.88 ± 0.01, (3 = 0.68 ± 0.01, z = 
1.65 ± 0.02, = 1.24 ± 0.01 and a = 0.42 ± 0.01 
(Fig. [5]). It may be noted that these independently 
estimated exponent values do not fit very well with 
the scaling relation j3 = v\\u. However, this type 
of scaling violation is also obseved in many active- 
absorbing transition cases [2fj| . 
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Model B 


Manna 


/3 


ID 
2D 
MF 


0.68 ± 0.01 
0.98 ± 0.02 


0.42 ± 0.01 
0.67 ± 0.02 
0.99 ± 0.01 


0.382 ± 0.019 
0.639 ± 0.009 
1 


z 


ID 
2D 
MF 


1.65 ± 0.02 
2.00 ± 0.04 


1.55 ± 0.02 
1.55 ± 0.02 
2.0 ± 0.04 


1.393 ± 0.037 
1.533 ± 0.024 
2 


a 


ID 
2D 
MF 


0.42 ± 0.01 
1.00 ± 0.01 


0.16 ± 0.01 
0.42 ± 0.01 
1.00 ± 0.01 


0.141 ± 0.024 
0.419± 0.015 
1 


v \\ 


ID 
2D 
MF 


1.24 ± 0.01 
1.00 ± 0.01 


1.90 ± 0.02 
1.20 ± 0.03 
1.00 ± 0.01 


1.876 ± 0.135 
1.225 ± 0.029 
1 



TABLE I. Comparison of critical exponents of this model with 
those of the conserved Manna model [27J • 



B. Model B 

1. Mean field case 

For the mean field case, we have studied for N = 10 6 , 
averaging over 10 3 initial condition. We numerically in- 
vestigate the phase diagram and the universality classes 
of the transition. In mean field case, the phase boundary 
seems to be linear starting g c — 1/2 for p — and ending 
at g c — 1 ior p — 1 (Fig. [3]), obeying g c — i(l+p). In this 
case, for p — 0, we find the critical point to be 0^= 1/2, 
and this is similar to the fixed energy sandpiles [25| . The 
critical exponents are the same along the phase boundary 
and they match with those of model A. 



IV. ANALYTICAL TREATMENT OF THE 
MODELS IN MEAN FIELD CASE 

In the mean field case the particles can jump from a 
site to any other, all with equal probability, or, equiv- 
alently, the underlying graph is fully connected. From 
the theory of zero range processes we know that the sta- 
tionary probability distribution is factorized in terms of 
single site quantities. If p % n is the probability that a site 
i has n particles, the average rate of outgoing particles 
is (u l ) — J2 n Pn v ( n ) n - By symmetry, we can work un- 
der the hypothesis that p l n — p n Vi, we have simply 
(it 1 ) = (u) = J2 n nv ( n )P n ' t na t a t stationarity is equal 
to the average number of incoming particles, since the 
density is fixed. A sequential dynamics would consist 
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FIG. 4. Simulation results for the case p = in 1-d, g c = 0.892 ± 0.001. (a) Variation of steady state density p a of active 
sites versus g — g c , fitting to ft = 0.42 ± 0.01. The inset shows the variation of p a with density g. (b) relaxation to absorbing 
state near critical point for different system sizes L, the inset showing the scaling collapse giving estimates of critical exponents 
a = 0.15 ± 0.01 and z — 1.40 ± 0.02. (c) Scaling collapse of p a {t). The inset shows the variation of p a (t) versus time t for 
different densities g. The estimated critical exponent is v\\ = 1.90 ± 0.02. The simulations are done for linear chains of size L 
(=N). 



of sequences of one-particle jump events and the master 
equation has the form: 

Pk (t) = -((«) + kv(k))pu(t) 
+6{k){u)pk-i( t ) + ( k + !M fc + 

(4) 

Multiplying by s k , summing over k we have the equation 
for the characteristic function G(s) — J2k Pk sk '- 

oo 

G = -(1 - s)((u)G(s) - J2 ™(»K* n-1 ) (5) 

n=l 

from which we have a self-consistent formula for the sta- 
tionary solution: 

-. oo 

G( S ) = —Y / nv(n)p nS n - 1 . (6) 
(11.) * — * 



Q_ 




FIG. 5. Behavior of order parameter for sequential and par- 
allel dynamics. For sequential dynamics, g c = 1 for both 
models. For the model with p = 0, the parallel dynamics 
gives g c = \. 



For the case (A) the solution of Eq. (j6]) is 

G(s) = 0(1 - g){\ -g + sg) + s6(g - lJe-fo-DC 1 -') (7) 
For the case (B) the solution of Eq.© is 

G(s) =0(l-g)(l-g + sg) + 0(g - 1) € -^± (g-x), (8) 

where x is the solution of g = ^ e _ 1 . In both cases below 
g c = 1 the activity is zero, i.e. the system falls into the 
absorbing state till this is present. The order parameter 
or the fraction of active sites as a function of g is depicted 
in Fig. [SJ On the other hand, the parallel dynamics of 
the model B with p = is particularly simple to analyze, 
since the number of particles on top of each site follows 
the discrete evolution equation: 



n t +i 



n l if n(t) y£ 1, 
n* + 1 if n(t) = 1. 



Then, we have the self consistent equation for the sta- 
tionary state: 

(5 -Pi)" , a , v (9 -Pi)"" 11 



(I -Pi) 



(n-1)! 



(9) 

From we which we can work out the characteristic func- 
tion 



G(s) = (1- Pl + s Pl )e- {9 - pMl ~ s) 



(10) 



From G'(0) = pi, we end up with an equation for x = 
g — pi, that, apart from the solution x = 0, has the form 
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x(l + xe x ) 
1 — e~ x + xe~' 



(11) 



Finally, the order parameter can be calculated consis- 
tently once we know G(s) from x (Eq. (jHJ), since p a — 
1 — G(0) — G'(0) and its behavior as function of g is de- 
picted in Fig. [5] 



Approximate analysis of the critical point and 
faster-is-slower effect 



To get insights on the value of the critical point as well 
as on the time to reach a frozen configuration below g = 
1, we will analyze the dynamical stability of the frozen 
phase probing it with a simple perturbation of the form: 



Pa = 1 - g + 5 
Pi = 9 - 2<5 
P2 = S 



(12) 
(13) 
(14) 



i.e. we pick a fraction 26 of particles and move them 
to already filled sites, neglecting the case in which some 
of them choose the same site. This perturbation evolves 
according to the equation |28[ 

5(t + 1) = S(t) -(l-g + S(t))(l - e - 2v{2Wt) ) 

+5{t)v^) 2 ^ 2vi2)m (15) 
whose solution scales, for large enough times, like 



5{t) oc 



2v{2) 1 - 



v(2) 



with a relaxation time 



log 1 - 2v(2)(l -^-g) 



(16) 



(17) 



that diverges in g c = 1 — -^p. This value coincides in 
a very good approximation with the numerical values of 
the critical point both for for the model A, g c — 3/4 and 
B, g c = ^(1 + p)- 

Moreover, for the model B the relaxation time at fixed 
density as a function of p 



-(P) 



1 



log [1 



2(1 -p) 



' l+P 
. 2 



9)} 



(18) 



has a minimum, optimal value r* = — i og [ g (2_ g )] m P* — 
<?, as it is shown in Fig. [6j This is a very simple example 
of the faster-is-slower effect, by which agents who choose 
a strategy which imply longer waiting times allow for a 
faster collective coordination. Note that waiting can be 
rational at the individual level, if an agent expects that 
others are leaving her site. 



B. Analysis of the finite size effects on the time to 
reach the absorbing state. 

Let's call ir the weight of the frozen configuration for 
given N sites and K = gN particles, its inverse can be 
an estimator of the time required to reach the absorbing 
state. We have: 




0.25 



0.5 
P 



0.75 



FIG. 6. Relaxation time r as a function of the waiting prob- 
ability p for the model B for a few values of g. It diverges at 
p = 1 and p = 2g — 1 and has an optimal minimum at p* = g, 
when r = r*. 
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N 

FIG. 7. Exponential dependency of timescale to reach the 
absorbing state, with the system size in mean field cases for 
Model A at g = 0.55 and Model B at g = 0.80. The data are 
averaged over 10 4 realizations. 



Since the process is zero-range, at stationarity P{n) — 
Elf Pi( n i)i then, if the graph is homogeneous we can 
leave the dependency on the site to get 



K N-K 

Pi Po 



Using Stirling approximation we have finally 



log7r 



-N 



g log A + (1-g) log (l—l 
Pi V Po 



(20) 



(21) 



i.e. it is exponentially decreasing with the system size. 
The Fig. [7] shows some typical numerical results for the 
average time to reach an absorbing configuration as a 
function of the system size. 



SUMMARY AND DISCUSSIONS 



7T= Yl P<jt) 

frozen 



(19) 



The dynamics of social systems is recently gaining more 
and more insights borrowing concepts and techniques 
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from statistical mechanics and by developing some toy 
models. The toy model paradigm [29[ helps us to under- 
stand complex emergent behavior of such socio-economic 
systems. In this paper we show how a crowd dynam- 
ics in a resources allocation game gives rise to a phase 
transition between an active and a frozen phase, as the 
density varies. In this respect, we have defined and stud- 
ied a class of models, where gN agents compete among 
themselves to get the best service from N restaurants of 
same rank, generalizing the 'Kolkata Paise Restaurant' 
problem. In the original problem, where density g — 1 
was far from its critical value g c , the relaxation time r, 
given by Eqn. (|3]), never showed any L = N x l d depen- 
dence. We recast these models in terms of zero-range 
interacting particles in order to get analytical insights on 
the systems' behavior. As long as g < 1, absorbing frozen 
configurations are present, and that can be reachable or 
not, depends on the underlying dynamics. We found the 
existence of a critical point g c above which the agents are 



unable to find frozen configurations. In the case in which 
the agents are moving if and only if they are unsatisfied 
(model B) with p = 0, they fail to reach satisfactory con- 
figurations if the density is above g c = 1/2. Strategies 
where agents wait longer (higher p) speed up the con- 
vergence, increasing g c and decreasing the time to reach 
saturation configurations (faster-is-slower effect). We in- 
vestigated numerically the phase transition in finite di- 
mensions finding a good agreement with the exponents 
of stochastic fixed-energy sandpile (Table. 01 p5l.l27l.[3Cl|. 
Thus, we have a simple model for resource allocation, 
which is solvable, and shows a variety of interesting fea- 
tures including phase transitions as in well known models. 
Further investigations are needed in order to understand 
the emergence of waiting strategies, i.e. how they could 
depend on agents' beliefs or learning processes, as well 
as it would interesting the study the case of heteroge- 
neous agents (with different strategies) and/or restau- 
rants (with different ranks). 
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